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We solve the cosmological evolution of warm dark matter (WDM) density fluctuations within 

■ the analytic framework of Volterra integral equations presented in the accompanying paper In 
the absence of neutrinos, the anisotropic stress vanishes and the Volterra-type equations reduce to 
a single integral equation. We solve numerically this single Volterra-type equation both for DM 
fermions decoupling at thermal equilibrium and DM sterile neutrinos decoupling out of thermal 

^ ' equilibrium. We give the exact analytic solution for the density fluctuations and gravitational 

potential at zero wavenumber. We compute the density contrast as a function of the scale factor a 
for a relevant range of wavenumbers k. At fixed a, the density contrast turns to grow with k for 
k < kc while it decreases for k > kc, where kc — 1.6/Mpc. The density contrast depends on k and a 
mainly through the product k a exhibiting a self-similar behavior. Our numerical density contrast 
for small k gently approaches our analytic solution for fc = 0. For fixed k < 1/(60 kpc), the density 
contrast generically grows with a while for k > 1/ (60 kpc) it exhibits oscillations starting in the 
radiation dominated (RD) era which become stronger as k grows. We compute the transfer function 

■ of the density contrast for thermal fermions and for sterile neutrinos decoupling out of equilibrium 
in two cases: the Dodelson-Widrow (DW) model and a model with sterile neutrinos produced by 
a scalar particle decay. The transfer function grows with k for small k and then decreases after 
reaching a maximum at fc = fee refiecting the time evolution of the density contrast. The integral 

^ ■ kernels in the Volterra equations are nonlocal in time and their falloff determine the memory of the 

i—^ ' past evolution since decoupling. We find that this falloff is faster when DM decouples at thermal 

[/3 , equilibrium than when it decouples out of thermal equilibrium. Although neutrinos and photons can 

■ be neglected in the matter dominated (MD) era, they contribute to the Volterra integral equation 
in the MD era through their memory from the RD era. 
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I. INTRODUCTION AND SUMMARY OF RESULTS 

In an accompanying paper [l| we provided a framework to study the complete cosmological evolution of dark matter 
(DM) density fluctuations for DM particles that decoupled being ultrarelativistic during the radiation dominated era 
which is the case of keV scale warm DM (WDM). In this paper, we solve the evolution of DM density fluctuations 
following the framework developed in ref. [l] . 

The new framework presented in ref. 1] and here is generic for any type of DM and applies in particular to cold 
DM (CDM) too. The collisionless and linearized Boltzmann-Vlasov equations (B-V) for WDM and neutrinos in the 
presence of photons and coupled to the linearized Einstein equations are studied in detail in the presence of anisotropic 
stress with the Newtonian potential generically different from the spatial curvature perturbations. 

In ref. [l| the full system of B-V equations for DM and neutrinos is recasted as a system of coupled Volterra 
integral equations. (Ref. [3l has recently considered this issue in absence of anisotropic stress). These Volterra-type 
equations are valid both in the radiation dominated (RD) and matter dominated (MD) eras during which the WDM 
particles are ultrarelativistic and then nonrelativistic. This generalizes the so-called Gilbert integral equation only 
valid for nonrelativistic particles in the MD era. 

We succeed to reduce the system of four Volterra integral equations for the density and anisotropic stress fluctuations 
of DM and neutrinos into a system of only two coupled Volterra equations. 

In summary, the pair of partial differential Boltzmann-Vlasov equations in seven variables for DM and for neutrinos 
become a system of four Volterra linear integral equations on the density fluctuations A(;m(?7, fc), A^(?7, fc) and 
anisotropic stress Yidm{^,k) ,T,^(ri,k) for DM and neutrinos, respectively. 

In addition, because we deal with linear fluctuations evolving on an homogeneous and isotropic cosmology, the 
Volterra kernel turns to be isotropic, independent of the k directions. As stated above, the k dependence factorizes 
out and we arrive to a final system of two Volterra integral equations in two variables: the modulus k and the time 
that we choose to be as 



y = a{r])/aeg^3200 a{f]). (1.1) 

We have thus considerably simplified the original problem: we reduce a pair of partial differential B-V equations on 
seven variables rj, q, x into a pair of Volterra integral equations on two variables: 77, k. 

The customary DM density contrast 5{-q, k) is connected with the density fluctuations Adm{v^ k) by 0] 



Arf„i(T?, fc) ^ 

Pdm [aeq + a{r])] ' 3200 



6{v,k)^ ^ Z^,, , a,,^^, (1.2) 



where pdm is the average DM density today. 

It is convenient to define dimensionless variables as 



Ho m \ aeq Vldr, 



where Ifs stands for the free-streaming length [1, [13) US] 5 ^^'^ comoving DM decoupling temperature and I*" is 
the dimensionless square velocity dispersion given by 

/"OO 

= / Q" /o'"(Q) dQ , while /o'"(Q) is normalized by I^"" = 1 . (1.3) 
Jo 
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Q is the dimensionless momentum Q = q/Td whose typical values are of order one. 

A relevant dimensionless rate emerges: the ratio between the DM particle mass m and the decoupling temperature 
at equilibration, 

_ maeq _ m f gd y 

= ^ - 4900 — j , 

gd being the effective number of UR degrees of freedom at the DM decoupling. Therefore, S,dm is a large number 
provided the DM is non-relativistic at equilibration. For ni in the keV scale we have £,dm ^ 5000. 

DM particles and the lightest neutrino become non-relativistic by a redshift 

ztrans + 1 = ^ ^ 1-57 X lo"^ ( T^) ' ^^r DM particles , z^,.^^^ = 34 for the lightest neutrino . 

Id KcV VlOO/ 0.05 eV 

(1-4) 

Ztrans denoting the transition redshift from ultrarelativistic regime to the nonrelativistic regime of the DM particles. 
The final pair of dimensionless Volterra integral equations take the form 

A{y,a) = C{y,a) + B^{y)^{y,a)+ C dy' [G^{y,y') ^{y\a) + Gl^{y,y') a{y' ,a)\ , (1.5) 



a[y,a)^C%y,a)+ dy' [r^{y,y') a[y' ,a) + I^{y,y') ^[y' ,a)\ , (1.6) 



with initial conditions A(0, a) — \ , (t(0, a) = | . This pair of Volterra equations is coupled with the linearized 
Einstein equations. 

The kernels and the inhomogeneous terms in eqs. (|1.5p - (|1.6p are given explicitly by eqs. (|2.12p - (|2.15p . (I2.16p - (|2.18p 
and (|2.7p - (|2.1ip . The arguments of these functions contain the dimensionless free-streaming distance l{y,Q), 



Jo 



The coupled Volterra integral equations (|1.5p - (ll.6|) are easily amenable to a numerical treatment. 

During the RD era the gravitational potential is dominated by the radiation fluctuations (photons and neutrinos) . 
The photons can be described in the hydrodynamical approximation (their anisotropic stress is negligible). The tight 
coupling of the photons to the electron/protons in the plasma suppresses before recombination all photon multipoles 
except 00 and Oi. (The Qi stem from the Legendre polynomial expansion of the photon temperature fluctuations 
e(r;,g,fc)i). 

00 and 01 obey the hydrodynamical equations [3] 

— + ^0i(,,«) = ^ , (1.8) 

^~^0o(^,5) = ^0(^,d) . (1.9) 

This is a good approximation for the purposes of following the DM evolution 

The photons gravitational potential is given in the RD and MD eras by (ref. |2'] and Appendix 



0(r;,fc) = 7/;(r;,fc) = 3^(O,fc)— ji -| , K^k^* , ,y* = / ^ = 143 Mpc , (1.10) 

2/ V V 3 / \ Qm Ho 

where ji(x) is the spherical Bessel function of order one. 

For redshift z < 30000 the kernel Ga{y,y') in eq. (|1.5p simplifies as 

M y,y'>0.1 Cdm K. y y' 



Gaiy,y' 



f^;- n [a {s{y) - siy'))] where, 

2/5 + y' 
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Il{x)^ Q dQ f^'^iQ) sin{Q x) and s(y) = -Arg Sinh 



( 



1 



) 



(1.11) 



In this regime z < 30000, y > 0.1, the anisotropic stress a{y,a) turns to be neghgible and eqs. (|1.5p - (|1.6p becomes a 
single Volterra integral equation. In the MD era this equation takes the form 



where the inhomogeneous term g(jj,a) contains the memory from the previous times y < 1 of the RD era. When DM 
is non-relativistic the memory from the regime where DM was uhrarelativistic turns out to fade out as 0.0002 
compared to the recent memory where DM is non-relativistic. 

The falloff of the kernel 11 [a (s — s')] determines the memory in the regime where DM is non-relativistic. We 
find that this falloff is faster when DM decouples at thermal equilibrium than when it decouples out of thermal 
equilibrium (see fig. [5]). This can be explained by the general mechanism of thermalization |25|: in the out of 
equilibrium situation the momentum cascade towards the ultraviolet is incomplete and there is larger occupation at 
low momenta and smaller occupation at large momenta than in the equilibrium distribution. Therefore, the out of 
equilibrium kernel n(a;) which is the Fourier transform eg. p. lip of the freezed out momentum distribution exhibits 
a longer tail than the equilibrium kernel. 

Neutrinos and photons can be neglected in the matter dominated era. However, they contribute to the Volterra 
integral equation in the MD era through the memory integrals over < < 1, namely the memory of the RD era. 

When the anisotropic stress (T{y, a) is negligible, eas. (|1.5p - (|1.6p reduce to a single Volterra integral equation for the 
DM density fluctuations Adm{y, a) when the anisotropic stress a{y, a) is negligible. We find the solution of this single 
Volterra equation for a broad range of wavenumbers 0.1/Mpc < fc < 1/5 kpc. 

At zero wavenumber k = the kernel of this Volterra equation vanishes and the DM fluctuations can be expressed 
explicitly in terms of the gravitational potential (/>. The gravitational potential at fc = follows solely from the 
hydrodynamic equations for the radiation combined with the regularity requirement at fc = of the first linearized 
Einstein equation. Namely, the gravitational potential </> is solely obtained from the radiation without specifying the 
sources of the DM and radiation fluctuations. Using this explicit and well known form of 0, (see e. g. ref. 0) the 
DM fluctuations are obtained at a = 0. The fact that the Einstein equations constrain their sources was first noticed 
in ref. [l6j in a completely different context. 

We depict in figs. [2] the normalized density contrast vs. y (the scale factor divided by aeq) for thermal fermions 
and sterile neutrinos in the Dodelson-Widrow (DW) model Q (both models yield identical density fluctuations for a 
given value of £dm)- Similar curves are obtained in the x niodel where sterile neutrinos are produced by the decay of 
a real scalar 1211 . 




At fixed y we find that the density contrast grows with fc for k < kc while it decreases for k > kc, where fcc — 1.6/Mpc. 
We find that the density contrast depends on a and y mainly through the product a y exhibiting a self-similar behavior. 
The density contrast curves computed numerically for small a gently approach in the upper fig. [2]our analytic solution 
for a = 0. For fixed a < 1, the density contrast generically grows with y while for a > 1 it exhibits oscillations starting 
in the RD era which become stronger as a grows (see fig. [2]). The density contrast becomes proportional to y (to the 
scale factor) at sufficiently late times. The larger is a, the later starts 5{y^a) to grow proportional to y (see fig. [2]). 
Also, the larger is a > 1, the later the oscillations remain. 

We depict in fig. [3] the transfer function for thermal fermions and sterile neutrinos in the DW model and for sterile 
neutrinos decoupling out of equilibrium in the x niodel. The transfer function grows with fc for small fc and then 
decreases after reaching a maximum at k ~ k^. 

We analyze in section [TTTl the system of two Volterra integral equations in the regimes where DM is in the transition 
from UR to NR and when DM is nonrelativistic. In sec. IIII CI we take the nonrelativistic limit of our system of 
Volterra integral equations in the MD era. This yields the Gilbert equation (plus extra terms). We find extra memory 
terms and different inhomogeneities arising from our system of Volterra equations. 

In section HVl we consider the zero anisotropic stress case where the system of Volterra equations reduces to a single 
Volterra equation. The numerical solution for the DM fluctuations in a broad range of wavenumbers is presented and 
discussed, as well as the transfer function and the analytic solution for zero wavenumber. 

We present in sec. |V]the distribution functions, main parameters and integral kernels for sterile neutrinos decoupling 
out of equilibrium and compare them to fermions decoupling with a Fermi-Dirac distribution. Finally, we present in 
sec. IVII the generalization of the Volterra integral equation for cold dark matter. 




2/> 1 



(1.12) 
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In the RD era where radiation fluctuations dominates the gravitational potential we derive in Appendix |X] a second 
order differential equation for the gravitational potential. We show that the solution of this differential equation is 
well approximated by the Bessel function of order one eq. ()1.10|) . 

We provide in Appendix [B] explicit and useful expressions for the free-streaming distance l{y,Q) [see eq. (|1.7p ] in 
the main relevant regimes. 

II. THE VOLTERRA INTEGRAL EQUATIONS AND RELEVANT PHYSICAL SCALES 

We recall here the pair of coupled Volterra integral equations derived in the accompanying paper [ij from the 
Boltzmann-Vlasov equations for DM and for neutrinos. 

A. Density fluctuations and anisotropic stress fluctuations 

In the companion paper [l[ we defined dimensionless density fluctuations Acim{y,Oi) and A^{y,a) and dimension- 
less anisotropic stress fluctuations a{y,a) factoring out the initial gravitational potential "0(0, A:) in order to obtain 
quantities independent of the k direction. These relevant quantities are expressed as 

Adm(y,a)= / - — £{y,Q) fo [Q) — tttt^n — - A^iv^o^)^ / -; — Q faiQ) 



= V(0,fc) , il^iy,!?) ^ tp{Q,k) '4;{y,a) and 1/1(0, a) = l, 

a{y,K) ^ ij{0,k) a{y,a) , a) = a) - -(/(y, a) , a-(0, a) = 0(0, a) - 1 . (2.1) 

We then introduce in ref. p| the combined density fluctuation A{y, a) 



1 A „^ , My) 



Adm{y, a) -I — A^(y, a) 



jdra 



, /e = ^ + i?.(0)~i?,(0) = 0.727 , A(0,a) = l, 

(2-2) 

where S^dm is the ratio between the DM particle mass m and the physical decoupling temperature at equilibration 
redshift Zeq + l^ 3200, 

= ^ = 4800 ^(^f- 5520 (i^) ' . ,2.3) 

We use here the dimensionless wavenumbers [l|, 



2 2 Td , , . . * _ / aeq 1 



K=kr]* and a = - — k — — 1^=^= k where r\* = J — — — 143 Mpc . (2.4) 

kdm Hq m y/tteq ^^dm V -"0 

Using A{y, a) and (f{y, a) in ref. [l| allowed to reduce the system of four Volterra integral equations into a the 
following pair of Volterra integral equations: 

A{y,a) = C{y,a) + B^{y) ^{y,a) + f dy' [Gc.{y.y') ^{y\a) + Gl{y,y') a{y' ,0)] , (2.5) 



a{y,a)=C%y,a)+ dy' [r^{y,y') a{y' ,a) + I^{y,y') ^{y' ,0)] , (2.6) 

JO 

with the initial conditions [l[ 

A(0,a) = l , (T(0,a) = ^ ~ ^ i?,(0) . 



We have in eqs. (l23)) - ((2:6)) 

C(.,a)==-^ 



a{y,Q) , Ry{y) 



+ {y,a) , C{y,a) = — + —— a {y,a) , (2.7) 

J3 J t,dm J3 
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Ga(y,y') 
Gl{y,y') 

iUy,y') 



-^NM) + ^N:^{y,y') 



7^N:{y,y')-^NZ{y,y') 



2/4 LCdm " 2/3 



^ C/„(jy,y') + ^t^r(2/,y') 



1 

^dm 



K{y,y')~^u:^{y^y') 



(2.8) 
(2.9) 
(2.10) 
(2.11) 



In eqs. (|2.5p - p.6p we can use ~ i?,/(0). The DM integral kernels and inhomogeneity functions in eqs. p.7p - ()2.1ip 
are given by 



<i{y,a) = / Q'^ dQ e(y,Q) 
Jo 

y £.dm bdm{y) 
No.{y,y')-- 



£(y,Q) 



/o"(Q)c-2™(Q) + 0(o)^ 

/o™(Q) [4Q' + 3(^,„y)2] , 



Jo 



^Ql{y,Q) 



dm 



Q2 dQe{y,Q)^ji [a Iqiy.y')] 



dQ 



e{y',Q) 



e{y',Q) 



Kiy^y') 

a"{y,a) = 

Uc.{y.y') 
K(.y,y') 



dm 



0" dQ J, [a Igiy^y')] e{y,Q) e{y\Q) 



00 



Q^ dQ 

n'^y'^Jo siy^Q) 



fr{Q)^cim{Q)+m 



dfS 



dm 



d\nQ 



J2 



^Qiiy,Q) 



(2.12) 
(2.13) 
(2.14) 
(2.15) 
(2.16) 



Q^dQ dfg 



dm 



^K^y'^Jo e{y,Q) dQ 



e{y',Q) 



Q' 

£{y',Q) 



{2 n[alQ{y,y')]-Z j3[alQ{y,y')]} , (2.17) 



00 ^4 



Q" dQ dfH 



dm 



bn'^y^Jo £{y,Q) dQ 



£{y', Q) {2 ji [a Igiy, y')] - 3 js [a Igiy, y')]} 



(2.18) 



The function (^„^{Q) determines the intial conditions. We have for thermal initial conditions (TIC) and for thermal 
initial conditions (TIC) Q 



1 din//?™ ( 1 dln/o" 

, — for thermal initial conditions (TIC) , - —r, — tt 

cUQ)-l 2 ^'-Q c-0(Q).J 2 ding 



for TIC , 



(2.19) 



2 for GIC 



-2 for Gilbert initial conditions (GIC) 
The neutrino integral kernels in eas. (|2.7l) - (|2.1ip and inhomogeneity functions are given by 

N:i^{y,y')^~8r,n[Kr{y,y')] , 



t/r(y,2/' 



24 /i' 



5 y"^ 



{2 ji [k r(y,y')] - 3 js [k r{y,y')]} 



(2.20) 
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a"^(y,a) = -21^ [l + 2 <^(0)] jo [k r{y,0)] 



a"'- '^(y, a) --6/3^ [1 + 2 0(0)] 



.12 [k r{y, 0)] 



y2 



(2.21) 
(2.22) 



The Volterra integral equations (j2.5p - (|2.6l) are coupled with the linearized Einstein equations derived in the accom- 
panying paper pj 



(f'iy,a) = [l+7^o(y)] cr(y,a)-7j-^ Admiy,a)-^^^^ A^{y,a)~2 R^{y) Qo{y,a) 

(2.23) 



2e 



dm 



Here, 



^ Pdmjy) 

~ Pr{y) 



Q 



and 



(2.24) 



jdm 



^[i + 0{e,^y')] , ^dmy<l, 



O 



£,dm y>5 



(2.25) 



B. Relevant scales in the ultra- relativistic and non-relativistic DM regimes 

The evolution of the DM fluctuations presented here is valid generically for DM particles that decouple at redshift 
Zd, being ultrarelativistic in the RD era and become non-relativistic in the same RD era. That is, the evolution 
presented here is valid as long as ^dm ^ 1 which is the case from eq. (|2.3p provided DM decouples ultrarelativistically 
deep enough in the RD era. 

The framework presented in this paper is general, valid for any DM particle, not necessarily in the keV scale. More 
precisely, the treatment presented here is valid for ^ 1 and: 

1»^!^ = ^= 3200 ^ 

Tdphys Td Zd Zd 

which imphes Zd ^ 3200 £^dm- The redshift at decouphng turns to be 

Tdphys_ ^ 15 Tdj>hys ( 9±\^ ^ (2.26) 

Td 100 GeV VIGO/ ^ ' 

where we used Td = {^/gdY'^ T„nb and r„„f, = 0.2348 meV. 

DM particles are ultra- relativistic (UR) for z > ztrans, ztrans being the redshift at the transition from ultra- 
relativistic to non-relativistic DM particles 

Ztrans + 1= — ^ 1.57 X 10^ ( . (2.27) 

Td keV VlOO/ ^ ' 

Then, they become non-relativistic (NR) for z < ztrans- In terms of the variable y [ea. (|l.ip ] the transition from UR 
to NR DM particles takes place around y ^ y trans while decoupling happens well before y trans by y ^ j/cj: 

ytrans = 1/^dm ^ 0.0002 , yd = 3200/zrf ~ 2 X 10-12 . 

Notice that modes that reenter the horizon by the UR-NR transition y ^ ytrans, have from eqs. (2.30) and (2.40) of 
the accompanying paper ref. [l|, wavenumbers 

dm 



idm _ 2 -y/ 1 



ytrans V* Ifs Ifs 

That is, when DM particles become nonrelativistic the free-streaming length Ifs is of the order of the comoving horizon 

El. 
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Universe Event 


redshift z 


_ a + 1 3200 

^ ae, 2+1 2+1 


DM decoupling 


rp 1 

-J - 1 fi in^^ '^^ { 3d\?. 

z-d r-u ru j 


yd =i 2 X 10-12 


neutrino decoupling 


^ 6 X 10" 


0.5 X 10-« 


DM particles transition from UR to NR 


^ 1-6 X 10^ ^ 


2/tr-an. = 0.0002 

10"^ <y< 0.01 


Transition from the RD to the MD era 


==: 3200 




The lightest neutrino becomes NR 


_ at; 

0.05 eV 


0.05 eV 

Vtrans 


Today 


20 = 


yo ~ 3200 



TABLE I: Main events in the DM, neutrinos and universe evolution. 



At decoupling, the covariant neutrino temperature, decoupling neutrino redshift and y variable are, 

Trf^ 0.17 10~^ eV , z2~6xlO^ and y"^ ~ 0.5 x IQ-^ . 
The lightest neutrinos become non-relativistic at a redshift 

-o^ and 7;^ -94 "'Q^ eV 

'^trans ~ ^"^ n nc U trans — '^^ 

0.05 eV mi/ 

Namely, neutrinos become non-relativistic in the MD era when their density as well as their fluctuations are negligible. 
Thus, we can treat the neutrinos as ultra-relativistic or neglect them. 



The neutrino and photon fractions of the energy density are defined in general as 

RAv) = —r-r = 7^, — rrrr- > ^-fi^) 



p{f]) n,. + a,(ri) VIm ' T - _|_ ^(77) VLm 



where Puij]), p-y^rj) and p{ri) stand for the neutrino, photon and total energy density, respectively. In the radiation 
dominated era fir ^ a{r/) and Ru{rj) + Ry{ri) — 1. The neutrino fraction changes after neutrino decoupling when 
the cosmic temperature crosses the e+ — e~ threshold, that is 01, 

0.727 , 4 X 10^ < z < 6 X 10^ 

RAv)^{OAl , 3200 <z<4xl09 . (2.28) 
, < z < 3200 

The quantity defined by eq. p.2p is dominated by the neutrino piece R,^{0) and takes the value 

~ i?„(0) = 0.727 . (2.29) 
In the MD dominated era both Ry{rj) and Rj{ri) become very small and can be neglected. 
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Range of Validity 




l{y,Q) 


UR DM particles 
^dm y < 1 

< y < 10"'^ 


Q 


idm y 1 


Transition regime from UR 
to NR DM particles 
10"'' <y< 0.01 




ArgSinh(«-^) 




NR DM particles 
0.01 <y< 3200 

Cdm y 3> 1 




-2.rgSm.(-L).log(i|.).i^4 




MD era 

y> 1 

NR DM particles 


^dm y 


' +logf«^'^-Ui « 



TABLE II: The different regimes ultra-relativistic (UR), transition and non-relativistic (NR) of the free-streaming distance 
Ky^Q)- Notice that the second (third) formula for l(y,Q) is also valid in the first (fourth) formula for < y < 10~^ (y 2> 1). 
In addition, the third formula of l{y, Q) for y matches for 2/^1 with the asymptotic behaviour of the second formula 
for l(y, Q). The precise behaviours of l(y, Q) are derived in Appendix iBl and given by eQS. (|B3|) - (|B8|) . When DM is UR l{y, Q) 
grows as the comoving horizon r]* y and thus free-streaming efficiently erases fiuctuations. When DM becomes NR l[y, Q) grows 
much slower and free-streaming is inefficient to erase fluctuations. 



We summarize in Table U the ranges of the redshift z and the variable y (the scale factor normalized to unity at 
equilibration) for the main events in the DM, neutrinos and the universe evolution. 

The free-streaming distance l{y, Q) is expressed by eq. (|1.7p . l{y, Q) can be expressed in general in terms of elliptic 
integrals. In the present case where ^dm ^ 5000 we find in appendix B excellent approximations to l{y,Q) in terms 
of simple elementary functions. We display the free-streaming distance l{y, Q) and the particle energy e{y, Q) for the 
different regimes in Table HT] It must be stressed that each of the four formulas displayed in Table [ll] match with its 
neighboring expression as discussed in appendix B. 

From eqs. (|2.3p and (|2.4p we obtain for the dimensionless variable a, 

a = 58.37 fl^y kpc . (2.30) 

m \ gd J 

In terms of a, the primordial gravitational potential eq.(3.17) in the accompanying paper [3] becomes, 

i/'(0,a) = — 3- — (kpc)-^G(a), (2.31) 



ai \aoJ \ m J gd 

where ao = 1.167 10"^ (keV/m) (100/.gd)^ and 

< G{d) G*{d') >=d{a-d') 



III. FROM THE ULTRARELATIVISTIC TO THE NON-RELATIVISTIC REGIME OF THE DM IN 

THE VOLTERRA EQUATIONS 



We investigate here the system of Volterra integral equations (|2.5|) - (l2.6p first in the transition regime for DM and 
then in the non-relativistic DM regime. 



10 



14 



l{y,Q'= 0.1)/Q vs. logi'o y 
l{y,Q = i)/Q vs. log^o i^- 
/(i/,Q = 10)/Q vs. logio 1/^- 



FIG. 1: The free-streaming length in dimensionless variables l{y,Q) divided by Q vs. logj^p y for Q — 0.1, 1 and 10. [Recall 
that Xfs = [t]* /£,dm) Q l{y,Q) [H]- We explicitly compute l{y,Q) in appendix 151 l{y,Q) is given in the different regimes by 
eqs.(|B2}, HB3[) . (|B8[) and (|B9|I . Notice that log^gy = corresponds to equilibration. We choose here £^dm ~ 5000. 



A. Transition Regime 



We consider here the coupled Volterra integral equations ()2.5l - (|2.6p ) in the transition regime from ultrarelativistic 
to non-relativistic DM particles 0.5 10~^ < y < 0.01 well inside the RD era where the neutrinos are ultrarelativistic 
and they have already decoupled. 

The second entry of Table |TT] the one-particle energy e{y, Q) ~ \J {^dm)'^ 2/^ + and the free-streaming length 
l{y, Q) applies now. Therefore, we have from eq. (jB3|) . 



i{y,Q) 



1-1 

16 VCdm 



Arg Sinh 



^dni y 
Q 



Q_ 



Q_ 



0{y^ 



i-A 

16 \Un 



Arg Sinh 



Cdm y 
Q 



These formulas are to be inserted in eas. (|2.12|) - (l2.18p for a{y,a), a'^{y,a), Na{y,y'), N^{y,y'), Ua{y,y') and 

u^{y,y')- 



B. Non-relativistic Regime 



We write here the Volterra integral equation in the non-relativistic regime 3200 > y > 0.01. 

The third entry of Table HIl for the one-particle energy e{y, Q) ~ y and the free-streaming length l{y, Q) applies 
in this case. Notice that the difference of the free-streaming lengths which appears in the integrand of the kernel 
Na{y,y') ea. (|2.14l) is now Q- independent because the DM particles are non-relativistic: 

^g(y,2/') = | [Ky,Q)~iiy\Q)]^Q [siy)-s{y')] , 
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where we used ea. (IB8|) . neglected terms O ([Q/^dm]^ logQ/^dm) in the free-streaming length l{y,Q) and 



s{y) = — Arg Sinh 



1 



1 



ds 



dy 2y VI + y 



(3.1) 



In the non-relativistic regime the kernels Na{y, y') and N^{y, y') in eqs. (|2.14p - (|2.15p both in the RD and the MD eras 
become, 

poo jpdm 

NM) = -K{y,y') = iUm)' yy' dQ ^ Ji {a Q [s{y) - s{y')]} . 

Jo dQ 

Integrating by parts df^"^ /dQ in the above integral leads to the simpler expression 

N^{y,y') = -Kiv.v') = -{Unf y y' Tl[a (s(y) ~ s{y'))] for y, y' > 0.01 , 

where 



poo 

U{x) = QdQ ft\Q) sin(g x) 
Jo 



(3.2) 



(3.3) 



That is, in the nonrelativistic regime the kernel Na(y, y') becomes the Fourier transform of the zeroth order momentum 
distribution f^"'{Q). Notice that 

n(0) = , n'(0) = 1 , 

where we used eas. (|1.3p and p.3p . 

In a similar way we obtain for the kernels Ua{y,y') and U^{y,y') [given by eqs. (l2.17p - (l2.18|) ] in the nonrelativistic 
regime, 



UM,y') = -Kiy,y') = - 



3 2/' 



5 K- y-^ Jo 



dm 



dQ 4- {2 ji [a Q [s{y) - s{y')]) - 3 js (« Q Hy) - s(y')])} • (3-4) 
dQ 



Upon integrating by parts this formula can be recasted in the simpler form 

3 a [s{y) - s{y')] y' 



U^{y,y') ^ ~U^{y,y') 



1^2 y3 



Q'dQC^Q) Jo (a Hy)-s{y')] Q) 



,h ja [sjy) - sjy')] Q) 



where we used also the angular integrals in the Appendix B of the accompanying paper ref. [l| . 

When the DM particles become nonrelativistic the anisotropic stress a{y, a) decreases fast as 1/(k y)^. For j/ a > 1, 
then 1/(k y)^ < 10~^ the anisotropic stress can be neglected. Therefore, for y a ^ 1 the system of Volterra equations 
3)) - (|2.6l) reduces to a single Volterra equation for the density fluctuations A(y,a). 



In this nonrelativistic regime where k y S> 1, e(y, Q) — £,dm y, the inhomogeneous pieces C(y, a) and (y, a) from 
eqs.jlHl), ((2l6t . ((23T|) . ((2:221) and ^J} become 



C(y,a) = --V / Q^dQ 



2 k 



foiQ) c'LiQ) + m 



dfi' 



d\nQ 



Jo 



My) 



[1 + 2 0(0)] jo[^r{y,0)] 



C'^(y,a) 



{k idmf y^ Jo 



Q^ dQ 



/o™(Q)4..(Q) + 0(o) 



dfS 



dm 



d\nQ 



J2 



-6i?,(y) [1 + 2,^(0)] 



J2 [k r{y, 0)] 



9 9 

y^ 



where we used ea.(|B8l 



(3.5) 



liy, Q) ~ ?^^^(y, Q) = 2 s{y) + log (8 Un/Q) 



(3.6) 
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The Volterra equation (|2.5p - (|2.6l) at y involves the integral over all y' in the interval < y' < y. Namely, we need the 
kernel Naiy,y') for all y' in < y' < y. Therefore, in the nonrelativistic regime y > yi — 0.01 we need the kernels 
with mixed arguments, where y' belongs to the transition or to the ultrarelativistic regime (0 < y' < j/i). 

We obtain from eq. (|2.14p for the mixed kernel 



N^y, y') = y dQ^^l^-Q 



log 



8 ^dn 

Q 



2 Arg Sinh 



Arg Sinh 



Crfm y' 
Q 



e{y',Q) 



e{y',Q) 



for 2/ > 2/1 = 0.01 , y' <yi = 0.01 . 



(3.7) 



The kernel Na{y,y') is proportional to (Cdm)^ when both y and y' are in the nonrelativistic regime [eg. (13.21) ] while 
it is proportional to (,dm when y is in the nonrelativistic regime and y' is in the transition or ultrarelativistic regimes 
[ea. (|3.7p ]. Namely, in the nonrelativistic regime, the memory of the transition regime and ultrarelativistic regime 
fades as 1/fdm 0.0002. 

In the MD dominated era y > 1, radiation (photons and neutrinos) can be neglected: Rv{y) — R-yiv) — 0. Once 
neutrinos are negligible, the anisotropic stress CT(y, k) becomes very small and can be neglected too. Therefore, we 
have for y > I dropping the neutrino contributions in eqs. (l2.2p . (|2.5I) and (l2.7l) - ([27 



Adm(2/, n) = -2^d„ A{y, k) , 



^dniiy, a) = a{y, a) + y ^dm bdm{y) (t>{y, a) 



MD era 

dy' 



(3.8) 



K 



\/i + y 



= N^{y,y') 4>{y',a) + K 



' dy' 



vr+^ 



jK{y,y') aiy',a) 



Notice that the integrals here cover the radiation dominated era < y < 1. That is, the memory of the neutrinos and 
photons during the RD era is preserved in the MD era. 

The linearized Einstein equation for the gravitational potential in the MD era (f)(jj,a) = 'il'{y,a) become from 



2/(i + y) ^ + i + y + ^(«y)' 



Hy^oi) 



1 



2^^ 



Adm{y,a) 



This equation can be solved as 



1 



2 idm y Ja 1 + 2/ 



dy' 



— — P4y,y') Adm{y',a) where I3t,{y,y') = — — - e 



kV3 



(3.9) 



(3.10) 



We compute in appendix A of the accompanying paper [l[ this integral in the asymptotic regime k 2/ 3> 1. We find 
at leading order from eq.(A2) of ref.[i|. 



Adm{y,a) , 



(3.11) 



2Cdm (k yY 

which corresponds to the Poisson's law. This result applies for k — ^dm a/2 ^ 1. 

The asymptotic expansion of the function bdmiy) for large y follows expanding the integral representation ea. (|2.13p 
in inverse powers of ^dm y- We obtain after calculation, 



bdmiy) 



rdm 



o 



{£.dm yr 



2 iUm y f 
where we used ea. (ll.3p . 

In the MD era we can approximate the inhomogeneous term a{y, a) in eq. p.Sp as 



(3.12) 



a(2/,a)=a*^^(2/,a) 



2^<im y 



QdQ 



a Jo I'^'iy^Q) 



d\nQ 



sm 



Qr^'iy.Q) 



(3.13) 



In eg. p. 81) we can approximate the kernel Na{y,y') for y' > 0.1 according to eg. p.2p and change the integration 
variable from y' to s' , defined as in eg. p.ip . 



-Arg Sinh -= 

Wy 



, y'^yis')^- 



ds' 



1 



sinh's' dy' 2 y' y/T+V 
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Then ea. (|3.8l) becomes 



Adm{y,o:) 6 

= g{y,a) + - 

y a 



>iv) 



(1) 



ds'n[a{s{y)-s')] Ad,n{y{s'),a) 



where n(x) is given by ea. p.3[) . we used eq. (|3.1ip and 

9{y, a) = + - / ^ 

y y Jo v^ + y 

The term proportional to a) in eq. p.Sp becomes neghgible using eqs. p.lip and p.l2p : 



= [N^iy, y') 4>{y', a) + {y, y') a{y\ a)] 



(3.14) 



(3.15) 



Urn bdmiy) 0(2/,") 



y>l 9 Adm{y,a) 



{k yf 



Notice that a^^{y, a) is exphcitly known. Once the Volterra equations (j2.5p - p.6p are solved from y — till y = 1 we 
expHcitly know g{y,a) from eg. p. 151) . Then, the Volterra equation p.l4p can be solved to find Adm{y,ct) for y > I. 



By setting 



Admiy,a) = y ddm{y,a) , 



(3.16) 



we obtain from the Volterra equation p.l4p the Gilbert-type equation for the density fluctuations valid when the DM 
particles are nonrelativistic. Notice that contrary to the original Gilbert equation only valid in the MD dominated 
era, our equation is valid for all y > 0.01, z < 300000 well inside the RD era, 



dd,niy,a) = giy,a) + - 



6 r^y^ ds' 



, ■ , 2^ n [a (s(2/) - s')] dd„i{y{s'),a) . 
OL Js{i) smh s 



(3.17) 



where g{y,a) is given by eas. (|3.13l) and (I3.15p . 

The memory piece in eq. (|3.15p turns to be smaller than the first term a^^^{y,a)/y by an order of magnitude or so 
as shown by numerical calculations. In addition, for y > 1 (MD era) we can neglect the logarithmic dependence in Q 
present in l^^{y, Q) ea. p.6p . We therefore have for a^^^{y, a) in the MD era from ea. p.l3p 



a''^{y,a) 



Urn y 



z{y) Jo 

where s{y) is defined by ea. (l3.ip and 



QdQ 



z{y) = a 



friQ)c'LiQ)+m 



dm 



d\nQ 



sin[z(y) 



For thermal and Gilbert initial conditions eqs. (|2.19p we can express a^^^{y, a) in terms of the kernel n(z) defined by 
ea. p.3p as 



y 



m + 2 



m<y)) 
z{y) 



n'(z(y)) 



+ i 



(3.18) 



where j — Q for thermal initial conditions and j — 1 for Gilbert initial conditions. 

Eqs. p.l4p and p.l7p have the same kernel as the Gilbert equation in the DM era pH . [l3 | as it must be. The 
inhomogeneous term g(y, a) differs to those in refs. [ill [l3|: this is so because g{y, a) takes into account the memory 
from the previous evolution of the fluctuations since the DM decoupling in the RD era considered here. In Appendix 
nil CI we derive a Gilbert- type equation in the MD era from eqs. (|3.13p . ()3.15p and (|3.17p valid in the MD era. The 
obtained Gilbert- type equation eq. p.lOp contains an inhomogeneous term corresponding to temperature perturbation 
initial conditions plus a memory term including the contributions from the RD era. 
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C. The Gilbert equation from the Volterra equation in the MD era. 

Eq. (|3.17p can be written as 

^ I a{v + vo) 2 ^ ^ "^^J 

n' [a {v + vo)] 



ddmiv, a)- - I dv' y{v') H [a (w - v')] ddmW , a) = idm \ [2 ^(0) + l] 



'9 



a 

U[a (v + vo)] 1 



a (v + Vo) 2 

where g = for thermal initial conditions (TIC) and g ~ I for Gilbert initial conditions (GIG), 



^ = s + l = l-ArgSinh( , log(8U0-l- 4.288... + -ln(^— j +-ln(^^j , (3.20) 



vi ~ 1 + s{l) and M[y, a] stands for the memory term containing the contribution of the gravitational potential and 
anisotropic stress from the RD era 

K dv' 

M[y,a] = - / [N^iy,y')^{y',a) + N:{y,y')aiy',a)] . 

V Jo VI + y 

On the other hand, the Gilbert equation in the MD era can be written as [HI, [3| 

SMD{y,a) / dv' yMDiv')'n[a{v{y) -v')] SMD[y{v'),a]= lL[av{y)] , (3-21) 







where L = G or T. In the notation of ref. (14] L = G or T corresponds to Gilbert or temperature perturbation initial 
conditions, respectively. 



T ( \ n(^) ... 1 

Ig{z) = , It[z) = - 



n'(.) + 2n(£) 



VMoiv) = ^ „ and VMoiv) = 1 ^ ■ (3.22) 

{l-vy yfy 

We see comparing eqs. p.l9p and p.2ip that the inhomogeneous terms are different. The inhomogeneities in the 
Gilbert equation (j3.2ip contain the functions Icict v) or It{<x v) while the inhomogeneities in the Volterra equation 
()3.19p contain Ixia {v + vq)] for TIC and a linear combination of Irla {v + vq)] and Ig[c( {v + vq)] for GIG. Namely, 
the argument v in the inhomogeneous terms containing the kernels 11 and 11' is shifted by the quantity vq given by 
ea. (l3.20p (A similar shift was noticed in ref. 18]). In addition, the inhomogeneous term M[y, a], memory of the RD 
era in the Volterra equation is necessarily absent in the Gilbert equation which only takes into account the MD era. 

In summary, choosing TIG at decoupling in the Volterra system of equations yields TIG a,t y — 1 for the Gilbert 
equation. On the contrary, choosing GIG at decoupling in the Volterra system of equations yields a linear combination 
of TIC and GIG at y = 1 for the Gilbert equation. One can thus say that TIC are stable under the evolution of the 
fluctuations. 

In order to complete the comparison of eq. ()3.19p in the late MD era with the Gilbert equation (|3.2ip . notice that 

— I ''^^ 1 and hence y(v) . 

We conclude that the Volterra integral equation (13.191) in the MD era is very close although not identical to the 
Gilbert equation in the MD era ea. (l3.2ip . The inhomogeneity in the Volterra integral equation (I3.19P for TIG has a 
factor in front and its argument is shifted by the constant with respect to the usual inhomogeneity in the Gilbert 
equation p.2ip . For GIG a linear combination of Gilbert and thermal initial conditions appear at ?/ = 1 for the 
Gilbert equation. In addition, the term M[?/, a] containing the memory from the RD era is present in the Volterra 
integral equation (|3.19p while such term is absent in the Gilbert equation p.2ip . 



1 = 1 - ArgSinh ( — ) 1 and hence y{v) = 3 as in eq.dSSS). 



IV. SOLVING THE VOLTERRA EQUATION FOR THE DM DENSITY FLUCTUATIONS (WITHOUT 

ANISOTROPIC STRESS) 

The formulation of the cosmological fluctuations evolution in terms of Volterra equations provides an efficient 
computational framework for both analytic and numerical treatment. In the following subsections we solve the 
Volterra equation for DM fluctuations in the absence of neutrinos, i. e. without anisotropic stress. First, we solve 
the Volterra equation numerically for a wide range of wavenumbers. Second, we find the analytic solution at zero 
wavenumber. 
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A. Numerical solution of the Volterra equation for a wide range of wavenumbers 



In the absence of anisotropic stress the radiation fluctuations during the RD era can be treated in the fluid approx- 
imation. Neglecting the DM gravitational potential in the RD era, the gravitational potential is given in the fluid 
approximation by (see ref. Q and Appendix 



%/3 



K y 



(l3{y,a) = (j){y,a) = ijj{y,a) ^ 3 ji — = , 0(0,0;) = 1, (4.1) 

ji(x) is the spherical Bessel function of order one. Notice that lim ji(x)/x = 1/3 . 

In the MD era we can neglect the gravitational potential produced by the radiation and take as gravitational 
potential the one sourced by the DM fluctuations eq. (|3.10p . In the absence of anisotropic stress, the DM fluctuations 
Adm(y,ct) from eq. (|3.8l) obeys the Volterra equation: 

fy dy' 

Adm{y,a) ^ aiy,a) +y ^dmbdmiy) <t)iy,a) + K / , , Na{y,y') (l>iy',a) , (4.2) 

Jo vi + y 

where a{y, a) is given by eq. (|2.12|) with ^(0) = 1. 

We present here the numerical solution of eq. (|4.2p where we smoothly match the gravitational potentials given 
by eas. (|4.1l) and (|3.10l) . The full numerical analysis of the system of Volterra equations (I2.5p - (|2.6p including the 
anisotropic stress will be the subject of future work where we will also compare our approach with the numerical 
solution of the ODE hierarchy of B-V equations 

For y > 0.01 the DM particles become nonrelativistic eq. (|4.2l) simplifies and takes the form of eq. (|3.17l) 

ddmiy,a) ^ h{y,a) + - [ . ^'^ n [a {s{y) ~ s')] dd„i{y{s'),a) , diy.a) = ^'^"(^'"^ ^ (4,3) 
a smh^ s' y 



where I ^{y, Q) is defined by eq. (|3.6 



hiy,a)^- ^ + - / -^N^iy,y')^iy',a) (4.4) 

y y Jo V^ + y 

and we have neglected the memory piece from the DM fluctuations in the UR regime but kept the gravitational 
potential of the photons which is dominant. Ea. (l4.3p is a closed integral equation of Volterra type that determines 
approximately d{y, a). We have checked numerically that eq. (|4.3p reproduces the solutions of the full Volterra equation 
()4.2|) within a few percent. 



MD 



{y,<^) , 1^ dy' 



From the numerical resolution of the Volterra equation (14.21) we find the normalized density contrast 



The density contrast S{y, a) is given by ea. (|1.2p and eq.(4.36) in ref. [l]. 

We depict in fig. [5] the logarithm of the absolute value of the normalized density contrast for fermions with 
S,dm = 5000 which corresponds to DM fermions in thermal equilibrium with m = 0.6736 keV and sterile neutrinos in 
the DW model with m ~ 1.685 keV (both models yield identical density fluctuations for a given value of ^dm)- In 
both cases we used thermal initial conditions. 

The density contrast generically grows with y for fixed a < 1 while it exhibits oscillations starting in the RD era 
for a > 1 which become stronger as a grows (see fig. As expected, the Jeans' unstability makes the density 
contrast proportional to y (to the scale factor) at sufficiently late times. The larger is a, the later starts S{y,a) to 
grow proportional to y (see fig. [2]). Also, the larger is a > 1, the later the oscillations remain. 

There exists a value a = ac — 0.1 determining the transition between two regimes. We separately display the plots 
corresponding to a < etc and a > ac- We find that for a < ac and fixed y, 6{y, a) increases for increasing a while 
the opposite happens for a > ac- Namely, S{y,a) at fixed y decreases for increasing a- We see from fig. [2] that the 
curves for logj^Q \6{y, a)\ vs. logj^Q y keep bending for decreasing a — >■ towards the a = curve. [The a = curve is 
obtained analytically in eas. (|4.8p and (|4.17p below]. 
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In figs. [2] we see that for both a < olc and a > ac varying a shifts the curves 5(j/,a) vs. y with respect to each 
other but keeping their form essentially unchanged. This property indicates that (5(?/, a) mainly depends on a and y 
through the product a y, namely in a selfsimilar manner. 

We have computed (5(j/,a) in the x-model for sterile neutrinos and found curves quite similar to the thermal case 



B. Analytic solution of the Volterra equation at zero wavenumber 

At a = (that is, fc = 0), the Volterra equation (|4.2|) (zero anisotropic stress) can be solved in close form since 
from eq. (|2.14p its kernel vanishes: Na=Q{y,y') — 0. 



The inhomogeneous term a{y,0) in the Volterra equation ()4.2|) becomes using eq. (|2.12p 

POO 

a{y, 0) - - y Cdm bdUy) + Q'dQ e{y, Q) ft\Q) c^dmiQ) , (4.6) 

Jo 

Therefore, the Volterra equation (14. 2p at a = simply relates the DM density fluctuations in terms of the gravitational 
potential as 



^dm{y, 0) = / g2 dQ e{y, Q) f^'^{Q) cL(Q) + V Um bd^iy) 
Jo 



0(2/,O)-l 



(4.7) 



More explicitly, for thermal (TIC) and Gilbert (GIG) initial conditions ea. (l2.19p . Adm{y,0) takes the form 



y Cdm bdrniy) 0(2/, 0) 



Adm(2/,0) 



for TIG , 



(4.8) 



for GIG 



We can obtain 0(?/, 0) solving the hydrodynamic equations for the radiation fluctuations p.8p - (ll.9|) together with the 
linearized Einstein equations for the gravitational potential in the /c — ^ limit 



-k^ <j)r{r], k) = 16 TT G a^ijf) p-ylrf) 



©r,o(?7, ^) + ^ Kl) 6r,i('7, k) 



(4.9) 



3ft(r,) -^ + fc2,/,(r,, fc)+ 3/1^(7,) i/-(r,,fc) = -4 7r Ga^ (77) 4 ^,(77) e,,o(r/, fc) + ^^^(77) 5d™(r/, fc) , (4.10) 
orj L J 

where (jiriv^k) stands for the radiation contribution to the gravitational potential and pdmiif) Sdmiv^k) for the DM 
fluctuations. Eq. (|4.10[) can be written in dimensionless variables as 

y [1+ Uoiy)] ^ + ^ (« y f 4>{y. a) + [1 + na{y)] ^{y, a) = -2 e,,o(y, a) - ^ TZoiy) hmiy, «) , (4.11) 
ay 3 2 



where TZo{y) is defined in eas. p.24p - (j2.25p and we used eq. (|A3l) 
Since the left hand side of eq. (|4.9l) vanishes at fc = we have 

k 



3/1(77) 



e,,o(?7,fc) + o(fc^) 



(4.12) 



We neglect radiation momenta higher than I = 1 thus neglecting the anisotropic stress and set ijj{ri, k) — (j){rj^ k). 
We have from ea. (|1.8p in the /c — > limit 



e,,o(??,O)-0(77,O) = c 

where c is a constant. 

The initial values ?/'(0, fc) and 6^,0 (0, fc) are related by the 77 — > hmit of ea. (|4.10p as 

7/.(0,fc) = -2e,,o(0,fc) , 



(4.13) 



(4.14) 



17 



log|Ad„i(y,0)| vs. logio V 
log |(5(?/, Q = 0.003)1 vs. logio y 
log |(5(?/, Q = 0.005)1 vs. logy, y 
log |(5(i/,Q = 0.01)1 vs. logic V 
log l(5(y, a = 0.03)1 vs. logi,, y 
\og\5{y,oi = 0.1)1 vs. logy, y 



12 
10 
8 
6 



logl5(y,a = 0.1)1 vs. log^, y 
\og\5{y,a = 0.5)1 vs. log^Q y 
\og\5{y,a = 0.7)1 vs. log^o y 
\og\5{y,a^ l.)l vs. login V 
\og\S{y,a = 3.)1 vs. log^, y 
\og\S(y,a = 10.)1 vs. logig y 




-1 







1 



FIG. 2; The ordinary logarithm of the normalized density contrast S{y,a) vs. log^Q y from the numerical resolution of the 
Volterra equation (|4.2p for DM fermions in thermal equilibrium with m = 0.6736 keV and for sterile neutrinos in the DW model 
with m — 1.685 keV. (Both models yield identical density fluctuations for a given value of ^dm). There is a value a — — 0.1 
determining the transition between two regimes. We separately display the plots corresponding to q < Qc and a > Oc- We see 
that for Q < Qc and fixed y, 5{y,a) increases for increasing a while the opposite happens for a > etc. S{y,a — 0) is plotted 
from the analytic solution ea. (|4.8|) - (|4.17|) for TIC. We see that the different curves have essentially the same shape and are 
shifted from each other in an almost selfsimilar manner indicating that 5{y, a) is mainly a function of k y. 5{y, a) generically 
grows with y for fixed a < 1 while it exhibits oscillations starting in the RD era for a > 1 which become stronger as a grows. 
5{y,a) becomes proportional to y at sufficiently late times. The larger is a, the later starts 5{y,a) to grow proportional to y 
and the later the oscillations remain. 
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up to small l/^(jm corrections. We thus find from eas. (|4.13p and (|4.14p . 



c=--0(O,O) = 3e,,o(O,O) 



(4.15) 



Inserting eqs. (|4.12p and (|4.13l) in eq. (|1.9p yields 



drj 



Qr.O 



+ 2e,,o(r/,0) = c, 



which in terms of the variable y becomes 



y 



dy 



y dTZoiy) 



2 [Uoiy) + 1] dy 



0r,o(y) = C . 



This first order differential equation can be resolved with the explicit solution 

2 

5 y^ 



.y3_y2^42/ + 8(l-v/yTT) 



(4.16) 



up to small corrections of the order 1/^dm because we set here TZoiv) = V [see eqs. (j2.24p - (|2.25p ]. 

Then, the gravitational potential follows from eas. (|4.13l) and (I4.15P and we recover the known expression 



0(2/, 0) = 



0r,o(y) 



1 



0(0,0) 2 2e,,o(0) 10 y3 
For zero or small redshift, eq. (|4.17p becomes 

y»i 9 , 1 



9 z/ + 2 - 8 + 16 (vyTT - 1 



0(y,o) 



10 5y 



(4.17) 



(4.18) 



It must be noticed that the known expression eq. (|4.17l) for the superhorizon gravitational potential (see for example 
ref. 0) follows here solely from the hydrodynamic equations for the radiation (|1.8p - (|1.9l) combined with eq. (|4.12l) . 
[Eq. (|4.12p follows from the first linearized Einstein equation ()4.9p in the fc — > limit]. Namely, 0(?/,O) and 6r,o(?/,0) 
are obtained without specifying the sources of the DM and radiation fluctuations. 

We can flnd the matter source of the superhorizon gravitational potential (/)(y,0) by inserting eq. ()4.17p in the left 
hand side of eq. (|4.10p for il}{-q, k) — (j){r], k) and k = Q. We obtain using the dimensionless variable y 



[l + 7^o(y)] 



dy 



0(y,o) = - 



2 + 27^o(2;)- 



y dTZoiy) 



dy 



er,o(y,0) , 



(4.19) 



Contrasting eq. (|4.19p with ea. (|4.1ip and using ea. (|2.25p implies a DM source 



Sdmiy,0)= 4- 



rfln TZq 
dhiy 



0r,o(y,O) 



4e^,o(2/,0) for Uny<l 



36^,0(2/, 0) for f dm y > 1 



(4.20) 



That is, combining the linearized Einstein equations with the hydrodynamic equations for the radiation requires for 
consistency a precise relation between the dark matter and radiation fluctuations. This is a consequence of the fact 
that the Einstein equations constrain their sources as was flrst noticed in ref. [l6| in a completely different context. 

Inserting ea. (l4.18p into ea. (l4.8p yields for the DM density fluctuations today 

1^ 



Adm{y,0) -idm y ' 



( 9 3 

5 5y 



1 + 



23 3 
I To ~ 5^ 



1 + 



for TIC 



for GIC 
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FIG. 3: The transfer function today T(7) vs. 7 = y Q for — 5000. The (red) solid line curve is for DM fermions in 
thermal equilibrium with m = 0.6736 keV and sterile neutrinos out of equilibrium in the DW model with m = 1.685 keV. The 
(blue) dotted line corresponds to sterile neutrinos out of equilibrium in the x-model with m = 0.7203 r~^/* keV where r is a 
coupling constant [see eq. (|5.1|l ]. That is, 0.9365 < m/keV < 1.665 [see eQ. (|5.4p ]. 7 is defined by eQ. (|4.23p . The presence here 
of a single maximum at 7 = 7c — 0.2 is consistent with the curves for S{y, a) and the value of Oc in figs. [2] Notice that the 
two transfer functions turn out to be very similar although they describe quite different dynamics. 



The normalized density contrast ea. (|4.5l) becomes today and for zero wavenumber, 



<5(2/,0) 



9 



i + o{- 

V 



for TIC 



10 



(4.21) 
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y 



1 + 1 - 

y 



for GIC . 

We depict log^g |5(y,0)| vs. log^Q y in fig. [5]for TIC. Fig. [5] exhibits this constant behaviour in 5{y,Q) for large y. 



C. The transfer function for the density contrast 



The transfer function at redshift z can be defined as the density contrast at redshift z > normalized by its initial 
value and then normalized by the whole expression at fc 0. That is, 



n.,«).f4 , r(,,o) = i. 

'5(2/, 0) 



The transfer function today becomes 



-1 ri rdra OH jdm 

T{a)^Yiu,T{y,a)^—^5{y,a) for TIC and T{a) ^ 5{y,a) for GIC , r(0) = 1 , (4.22) 
and we used eq. (|4.2ip . 
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We plot in fig. [3] the zero redshift transfer function T(7) vs. 7 for ^dm = 5000 and TIC. The (red) sohd line curve 
is for DM fermions decoupling in thermal equilibrium with m — 0.6736 keV and sterile neutrinos out of equilibrium 
in the DW model with m — 1.685 keV. The (blue) dotted line corresponds to sterile neutrinos out of equilibrium in 
the x-model with m ~ 0.7203 r^^/** keV. That is, 0.9365 < m/keV < 1.665 [see ea. (|5.4p ]. The variable 7 is defined as 




(4.23) 



We find that the transfer functions have a single maximum at ac consistent with the behaviour of h{y^a) in figs. [2l 

Notice that TW) grows fast with 7 for < 7 < 7c as we see from fig. [3l Previous calculations of the transfer 
function in refs. |12l.ll3| and with better precision in ref. (l3 | only exhibit the portion of T(7) where it decreases with 
7- 

The transfer function computed solely in the MD era monotonically decreases with 7 for growing 7 > 0, • The 
new piece of T(7) increasing with 7 for < 7 < 7c comes from the behaviour of the DM fluctuations in the RD era 
computed here, ac and 7c correspond here to a wavenumber fee ~ 1.6/Mpc. 

The transfer function for DM fermions decoupling in thermal equilibrium and for sterile neutrinos out of equilibrium 
in the x-model turn out to be very similar as seen from fig. |31 



V. FERMIONS IN THERMAL EQUILIBRIUM AND STERILE NEUTRINOS OUT OF EQUILIBRIUM 

The sterile neutrino is a serious candidate for WDM 0, [2l| . The freezed out distribution of sterile neutrinos turns 
to be out of thermal equilibrium in most models (22j . We consider here two sterile neutrino models for illustration. 
The Dodelson-Widrow model (DW) and the x model of ref. [2l| . 

The freezed-out DM distributions are given by 

DW model: /^^'^(Q) = \ /q - 0.043 keV , x model : /o'^(Q) = r /o^"(g) , 0.035 < r < 0.35 , 

(5.1) 

where t is a coupling constant and the normalized DM distribution function for the x model [9|, |23| takes the form 

/o*"(Q) = ■ (5.2) 

3 c(5) 0rg n2 

The normalized DM distribution in the DW model ^ is identical to the normalized Fermi-Dirac distribution 

/o*"(Q) (5.3) 

The simple formula eqs. ()5.1|) - (|5.3|) for the DW freezed-out distributions were given in and are widely used in the 
literature. A more sophisticated freezed-out distribution is derived in ref. (2^ . 

We plot in fig. S] the normalized distribution functions /□'"(Q) for fermions in thermal equilibrium (which is 
identical to the DW model) and for sterile neutrinos out of thermal equilibrium in the x model. 

We find from eqs. (|2.3|) and ()5.ip the values for the parameters and Ndm in the three DM fermion models 
considered here: 

Cfi' = 6721 (.gdm)^ (5!^)' , CdT = 2355 (5dm)^5^ , = 6146 (ff.„ r)^^^) ^ , 

N^^ = 1.805 , ATj^J^ = 0.07765 — , iV,^„ = 1.380 r . (5.4) 

m 

Notice that the out of thermal equilibrium distribution is larger than the equilibrium distribution for small momenta 
Q < 2 while the opposite happens for Q ^ 2. This can be explained by the general mechanism of thermalization: 
the momentum cascade towards the ultraviolet j25j . The distributions out of equilibrium therefore display larger 
occupation at low momenta and smaller occupation at large momenta than the equilibrium distribution. 
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Thermal FD and DW 


X model 







TABLE III: The normalized momenta J^" defined by eq-fOJ. Notice that = 1. 



^ Fermi-Dirac distribution 

sterile neutrinos in the x model 



FIG. 4: The ordinary logarithm of the normalized distribution functions /o™((5) vs. Q for fermions at thermal equilibrium 
(which is identical to the out of thermal equilibrium DW model) and for sterile neutrinos out of thermal equilibrium in the x 
model. 



We display in Table IIIII the momenta defined by eq. (|1.3p for the normalized DM distributions considered in 
this paper. 

For fermions decoupling ultrarelativistically at thermal equilibrium (and in the DW model of sterile neutrinos) the 
normalized freezed out distribution function is given by eq. (|5.3p . and the kernel n(a:) for the non-relativistic regime 
eg. p.Sp can be expressed as 

3C(3)^^(n2+a:2)2 " 

This kernel decreases for large argument x as 

For sterile neutrinos out of thermal equilibrium in the x model the kernel Tl{x) for the non-relativistic regime ea. p.3 
can be expressed as 



,/9 \/ (n^ + x^)2 + 3 n — 

= — iT^ ■ (5-5) 

3 C(5) ^1 n2 {rfi + x^)^ 
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This kernel decreases for large argument x as 



y2C(5/2) 1 



,3/2 



p5/2 



3C(5) 

We plot n^^(a;) and 11'^ (x) as functions of x in fig. [S] 

W^{x) has a longer tail than n-'^^(a;) due to the higher occupancy of the low Q modes in the out of equilibrium 
momentum distribution. The out of equilibrium kernel 11-^ [x) therefore exhibits a longer memory than {x) . 
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FIG. 5: The kernel n(a;) defined by ea. H3.3|) as a function of x for fermions in thermal equilibrium (which is identical to the 
DW model) and for sterile neutrinos out of thermal equilibrium in the x model. 



VI. VOLTERRA INTEGRAL EQUATIONS FOR COLD DARK MATTER 

All the framework of this paper easily generalizes to cold dark matter (CDM), that is DM particles with mass 
beyond one GeV that decouples nonrelativistically. For CDM the parameter S^dm is much larger than for WDM. 
Typically, 

_ m aeq _ 11 m 5 GeV 
- -jr- - 100 GeV Td,phys 

where we take 5 GeV as reference value for the physical decoupling temperature Td^phys of CDM. Other values for 
Td,phys appears in the literature according to the particle physics model chosen but ^cdm turns out to be very large 
in all cases (and much larger than for WDM where £,wdm ~ 5000). 

CDM decouples being (by definition) nonrelativistic thus we have at decoupling e(jjd,Q) = Ccdm Vd 2> 1. In 
addition, we have for the ratio of cdm to radiation densities TZoiu) = 2/ at all times after decoupling. 

CDM decouples at thermal equilibrium with a normalized Boltzmann distribution function [lo| 

/o''"(Q) = \/ where x = ^cd™ J/rf » 1 • (6.1) 

For DM particles decoupling nonrelativistic we have instead of eq. (|2.ip for DM decoupling ultrarelativistic, 

J 47r tp{yd,K,) 
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and therefore at the initial (decouphng) time y = yd 

/•OO 

Jo 

where we used the initial conditions for the BV distribution function discussed in ref.[l| 

'i'dm{yd,Q,K) = ■4'{yd,K) c°„(g) . 
From eqs.dHH) and we find as initial CDM fluctuations for TIC 

Acdm(2/d, = 2 j ^ ^2 ^ " 

For CDM the free-streaming distance is defined similarly to ea. (|I.7p as 
ry dy' _ dy' 

Jy^ y' \/i + y 



(6.3) 



= = lcdm{y) , thus, lcdm{y) = In — + 2 s(y) 

' ^ yd J 



(6.4) 



The free-streaming distance turns to be independent of Q as it should be because CDM particles are very slow. 

Let us consider zero anisotropic stress here, for simplicity. Thus, the evolution of the CDM fluctuations is given by 
the single Volterra equation (j4.2p . Since CDM particles are always nonrelativistic with ^cdm y ^ 1 the coefficients 
and kernel in ea. (|4.2l) take the form 



acdm{y, a) 



2 ^cdm y 

al{y) Jo 



QdQ 



dm 



d\nQ 







sin 





bcdmiy) = 3 , N^{y, y') = -{U,n? 2/ y' H [a {s{y) - s{y'))] , 



(6.5) 



where we used eqs. (|2.12p . (|2.13p . and (|3.2p . The Volterra integral equation for CDM takes thus a form similar to 
eg. (14:2)) for WDM 



(6.6) 



fy dy' 

Acdm(y,a) = acrfm(y,a) + 3 Cdm y (^(y,") + « / , , Na{y,y') (t>{y',a) and also 
Acdm{y, a) = Qcdmiy, a) +3 ^dm y 0(y, a) - 2 (^cdm)^ ny I . , 4 , n [a (s(y) - s')] (t>{y' 



sinh s' 

At y = yd, Icdvaiyd) = and eqs. ilOl) and (|6?5|) yield 

acdmiyd, a) = Acdm{yd, k) - 3 X . 

Therefore, the Volterra equation (|6.6p is identically satisfied at y ^ yd since (f)(yd,a) ~ 1. 

The whole section ITV Bl translates to the CDM case. Since for CDM ^cdm y ^ 1 eq. (|4.20p results for all y 

Scdrn{y,o) = 3 6^,0 (y,o) , 

as it must be. 
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Appendix A: The gravitational potential in the RD era 

During the RD era the gravitational potential 0(77, a) is dominated by the radiation (photons and neutrino) fluctu- 
ations. Neglecting the anisotropic stress, the following equations relate the gravitational potential with the first two 
radiation momenta 0, Q 



—k^ (j){ri, a) — IGn G a^{ri) Priv) 

— h k 6^,1(77,0;) = — , 

ar] drj 



(Al) 



Here 0r,o('7;Ck) and Or,i{r],d) are the first two momenta of the radiation temperature field 

0r,o('7,a) = Rjiv) 0o(fy,") + R-y{v) No{ri,a) and iVo(?7, a) = A,,(y, k) 0(0, a) . 

The infinite hierarchy of equations arising from the Boltzmann-Vlasov equation for radiation and matter, has been 
truncated to the first two equations 2, 3]. h(ri) = dlna/drj stands for the Hubble parameter and Prirj) — fir Pc/o,*{r]) 
for the radiation density. 

Eliminating 0,.^o(^jcf) among eqs. (|JT|) yields 



(j){r],d) 



1 



fc2 



a^('7) Pr(?7) 
0(?7,a) = . 



= 



16 TT G a'^{rj) pr{i]) 

It is convenient to use the variable y defined in ea. (|l.ip instead of the conformal time 77. We find in terms of y 

3 1 



a iv) Priv) 



I vr G 77*^ 



(A2) 



(A3) 



and eqs. (jA2[) read 
dQr 1 



1 



dy 1 + y \y 



1 1 K 



dQr,! 1 ^ . ^x 

— h - B,..i(y,a) = 

dy y 

Eliminating now 0^,1 (77, cf) we have 



^VT+y 



9 2 



3 Vl + y 
(j>{y,d 



9 9 

n y 



dy 



2 2 



(t){y,d) 



(A4) 



z 1 + 72^ 



y^ 



9 ? 

K y 

6 



dy 



(t>{y,d) = 



(A5) 



Taking the y derivative of this equation and replacing dQr,i{y^d)/dy and Qr,i{y,d) from eqs. (IA4[) and (IA5p respec- 
tively, we get the second order differential equation for the gravitational potential ^(77, d): 



-T^ + - R^iy) -1- + ^ ^>^iy) 't^^y^ ") = . 

dy y dy y^ 



(A6) 



where, 



R.iy) = 1 



4(1-Hy) l-(-K2y2/g l + |y + ^2y2/g 



1 + 1 (l-«^yV6)-«SV36 
(1 + y) (l + A«2yV6) 



(A7) 
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FIG. 6: The gravitational potential (j){y,a) and (jPiy^a) vs. logjQ for a — 0.1, 1 and 10 defined by eqs. (|A8[) and (|A9|) . We 
see that (^''(C) is a very good approximation to <f>{y,a) in the whole range of ^. 



In the radiation dominated era and for Ky ^ 1 this equation reduces to 



(f(j) AdS 1 9 , , 

dj/"^ y dy -J 



whose solution in terms of Bessel functions is given by eq. (|4.1 



We solve numerically ea. (|A6p both in the radiation dominated and in the matter dominated eras. We plot in fig. 
inithe normalized gravitational potential 

<^(y,a)^-^^^ , ^(0,a) = l, (A8) 
as a function of logj^Q C for a = 0.1, 1 and 10 as well as the function 

/(C) , C-^ , 0°(O) = 1. (A9) 

We see from fig. [5] that the function 4>^{Q is a very good approximation to 0(y, a). 

Appendix B: The free-streaming length in the different regimes. 

We evaluate now the integral eq. (|1.7p 

i{y. Q) = r ' (Bl) 

[l + y') \y'^ + {Q/Unf 
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in the different regimes depicted in Table III. This is an elhptic integral that can be expressed in terms of the standard 
incomplete elliptic integrals of first kind jlQl] 



1 + 



where 



2p^ = 1 



[F{v{0),p)-F{^{y),p)] 



1 + [Ql^d^nf 



(B2) 



Taking into account that ^rf,„ 5000 ^ 1 we can express this integral quite acurately in terms of elementary functions. 
• For y < 0.01 ^ 1 we can expand 1/^/T+l/ in ea. (|Bip in powers of y' and obtain 



V 



dy' 



1-A Q_ 

16 VCrfm 



+ i^dmf y 

1 

Arg Sinh 



^-\y'^\v'^ + o{y'^) 



id-m y\ 1 

' 2 



8 



Q V Q 



+ Oiy'). (B3) 



Notice that in this range of y and for typical Q ^ 1 the arguments in eq. (|B3l) can go from Q ^ £^dm y till 
Q ^ Cdm y- In the case ^dm y ^ Q £.dm this formula simplifies as 



i{y,Q) 



y 

Q 



• For y > 0.01 it is convenient to split the integral eq. (jBl|) in two pieces: 

dy' 



where 



;(0O,Q) = ^dri 



y Vi^ + y'W + iU^)'y"] ' 
dy 



(B4) 



In order to obtain the asymptotic expansion of l{oo, Q) for Q/^dm — > it is convenient to change the integration 
variable as 



y=t-l- 



1 



2 £,dm ) t~\ 



and l{oo,Q) becomes 



l{oo, Q) = 



dt 



1 



1+: 



'1 



iQ/S,dm) 
4t (t- 1) 



(B5) 



Expanding the integrand of ea. (jB5[) in powers of {Q /(,dmY ^^d integrating term by term yields the asymptotic 
expansion 



?(oo, Q) 



16 Ve* 



loa 



^dm\ 1 Q 7 / Q 



Q J 2 Cd,n 16 \^d 



O 



(B6) 



The integral in the second term of eq. (jB4[) is evaluated by expanding the integrand in powers of {Q/S,dm) 



idn 



dy' 



dy' I f Q 



^[i + y')[Q^ + {UnYy'^] Jy y'VT+V 2 \UnJ Jy y'-'VT+V 



dy' 



O 



Q 
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1 



16 ve 



Q 



Arg Sinh 



1 



Q 



O 



Q 



Collecting together eqs.(jB4]), dBB]) and dBT]) yields 

2" 



Ky,Q) 







1 - 


- (■ 


16 V 


7 




16 





Q_ 

ld7n 
2 



loe 



- 2 Arg Sinh 



1 



1 



8 \y £,drn 



(B7) 



(B8) 



• For y >• 1, ea. (|B8|) becomes 

Z(y,g)~log 







2 ^dm 



■^dm 



log 



We have thus derived the four entries for l{y, Q) in Table III. 

Moreover, ref. [l^l provides asymptotic expressions for the incomplete elliptic integral of first kind F{ip{y), k) valid 
for k'^ <C 1 (therefore Q Cdm) which are uniform in ip. The first order asymptotic expression takes the form j20| 



1 / 1 
1 + - (1-- 

p 



log 



-+\/- 

p \ u 



log 



1 + {Q/UmY 

2^ 2\u-p\ 



u = 1 + y . 



(B9) 



We display in fig. [l]Z(y, Q) for Q = 0.1, 1 and 10. Large values of Q get suppressed in the integrals by the decrease 
of the distribution function /o™(Q). 
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